function heat_step_movie

%  plots exact solution of the heat equation
%       D*diff(u,x,x) = diff(u,t)   for xL < x < xR, 0 < t < tmax
%  where
%      u = 0  at x=xL,xR  and  u = step(a,b)  at t = 0

co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

frames=4;

v = VideoWriter('/Users/mark/Desktop/heat_step_movie','MPEG-4');
open(v)

% get(gcf)
set(gcf,'Position', [1 996 692 349])
subaxis(1,1,1,1,'MT',0.08,'MB',0.1,'MR',-0.01,'ML',0.02,'P',0.04)


% set parameters
N=80;
M=80;
%M=10;
tmax=0.5;
xL=0;
xR=1;
a=1/3;
b=2/3;

% generate the points along the x-axis, x(1)=xL and x(N)=xR
x=linspace(xL,xR,N);

% generate the points along the t-axis, t(1)=0 and t(M)=tmax
%t=linspace(0,tmax,M);
t(1)=0;
dt=linspace(-4,0,M);
for it=2:M
    t(it)=10^dt(it);
end

nmodes=100;
U=zeros(N,M);
for j=1:M
    for i=1:N
        s=0;
        for ii=1:nmodes
            an=2*(cos(a*pi*ii)-cos(b*pi*ii))/(pi*ii);
            s=s+an*exp(-pi*pi*ii*ii*t(j))*sin(pi*ii*x(i));
        end
        U(i,j)=s;
    end
end

for j=1:M
    clf
    subaxis(1,1,1,1,'MT',0.08,'MB',0.1,'MR',-0.01,'ML',0.02,'P',0.04)
    axis([0 1 0 1.1])
    hold on
    
    say=['Solution of Diffusion Equation:  $u(x,0)$ shown below'];
    title(say,'FontSize',14,'FontWeight','bold','Interpreter','Latex')
    
    xlabel('x-axis')
    ylabel('Solution')
    
    box on
    
    if j==1
        plot([0 a],[0 0],'b','LineWidth',1.6)
        
        legend({' u(x,0)'},'Location','NorthEast','AutoUpdate','off','FontSize',14,'FontWeight','bold')
        
        plot([a a],[0 1],'--b','LineWidth',1.6)
        plot([a b],[1 1],'b','LineWidth',1.6)
        plot([b b],[0 1],'--b','LineWidth',1.6)
        plot([b 1],[0 0],'b','LineWidth',1.6)
    else
        plot([0 a],[0 0],'--b','LineWidth',1)
        plot(x,U(:,j),'-b','LineWidth',1.6)
        
        legend({' u(x,0)',' u(x,t)'},'Location','NorthEast','AutoUpdate','off','FontSize',14,'FontWeight','bold')
        
        plot([a a],[0 1],'--b','LineWidth',1)
        plot([a b],[1 1],'--b','LineWidth',1)
        plot([b b],[0 1],'--b','LineWidth',1)
        plot([b 1],[0 0],'--b','LineWidth',1)
        
    end
    
    yticks([0 1])
    yticklabels({'0','1'})
    
    xticks([0 1/3 2/3 1])
    xticklabels({'0','1/3','2/3','1'})
    
    say=['t = ',num2str(t(j),'%4.3f')];
    text(0.85,0.7,say,'FontSize',20,'FontWeight','bold')
    
    set(gca,'FontSize',15,'FontWeight','bold')
    
    % make movie frame
    F = getframe(gcf);
    for i=1:frames
        writeVideo(v,F);
    end
    
    hold off
    
end






